/******************************************************************
CLASS: lattice_1DTRI
DESCRIPTION:
KMC simulator with 1D triangular lattice considering three phenomena:
  1. diffusion
  2. confliction detection
  3. deposition

INPUT:
  1. W: deposition rate, layer/s
  2. T: temperature, K

OUTPUT:
  1. alpha_n and beta_n, n = 0,1,2,...,MODE:
  2. Rough2: roughness square, layer^2
  3. SOR: site occupancy ratio
  4. h: thickness, layer
  5. slope: root mean square slope

AUTHOR: Xinyu Zhang based on Gangshi Hu's code

CREATED: 4/22/2009

UPDATES:
7/29/2009
1. Introduce T as a control variable
2. Allow zero initial deposition rate
3. Use an algorithm in update() that can precisely indicate the deposition duration
4. Terminologies update: dens->SOR, v_->rec_, ht->H,
5. Intrudec structure lattice_state for the interface of lattice output, which is also used by the controller
*******************************************************************/

#ifndef LATTICE_1DTRI_H
#define LATTICE_1DTRI_H

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
using namespace std;


// #define DIFF_E0				(0.9) // For multiscale simulation
#define DIFF_E0             (0.6)
#define DIFF_R0             (1e13)
#define INT_COL             (PI/LatticeSize)	// interval between two (semi-)columns
#define M_INI				(100)					// initial lattice layer number
#define M_INC				(100)				// increasement of lattice layer number
#define SURF				(100)				// surface layer numbers
#define RELAX				(1)
#define MODE				(100)

#define EV2KB				(11604.505)			// the value of ev/k_B, the unit is K (Temperature)
#define H_ML				(0.8660254037844386)//sqrt(3)/2
#define PI					(3.141592653589793)
#define RECP_SQRT_PI		(0.56418958354776)	// 1/sqrt(pi)
/******************************************************************
Definition of bits:
the first three bytes are undefined and the default values are 0x00
the last byte is defined according to the occpancy and the neighbor state of the site:
0th bit (lowest): occupancy of the site
1st-6th bits: occupancies of the neighbors
	1-left; 2-bottom left; 3-bottom right; 4-right; 5-top right; 6-top left;
7th bit (highest): not defined
Occupancy: 1-occupied; 0-unoccupied.
*******************************************************************/
#define ST_INT					(0x0000007f)	// interior site
#define ST_FIX					(0x0000001f)	// fixed layer site
#define GET_NB(lyr,st)			((pp_lattice[lyr][st]&MASK_NB)?0x1:0x0)
#define GET_NB_NUM(lyr,st)		(NB_NUM[(pp_lattice[lyr][st]&MASK_NB)>>1])

// Example: GET_SITE(lyr,st,0) returns 1 if the site is occupied, 0 if unoccupied
#define GET_SITE(lyr,st,n)		((pp_lattice[lyr][st]&MASK_ST[n])?0x1:0x0)

#define SET_SITE(lyr,st,n)		(pp_lattice[lyr][st]|=MASK_ST[n])
#define RESET_SITE(lyr,st,n)	(pp_lattice[lyr][st]&=~MASK_ST[n])
const unsigned int	MASK_ST[7]=	{0x1,0x2,0x4,0x8,0x10,0x20,0x40};
const unsigned int	MASK_NB  =	0x7e;
const unsigned int	NB_NUM[64]=	{0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6};
// *******************************************************************
struct lattice_state
{
	double	SOR;
	double	pSOR;
	double	H;
	double	R2;
	double  alphas[MODE];
	double	betas[MODE];
	double	w;
	double  top;
	double  M2;
};

struct OneRecord
{
	double time;
	double T;
	double w;
	double H_mean;
	double H_var;
	double R2_mean;
	double R2_var;
	double SOR_mean;
	double SOR_var;
	double pSOR_mean;
	double pSOR_var;
	double M2_mean;
	double M2_var;
};

struct KMC_settings
{
	double t_start;
	double t_final;
	double t_int;
	int    lyr_ini;
	int    lyr_inc;
	int    lyr_surf;
};

struct KMC_model
{
	double diff_eng;
	double diff_R0;
	double mode;
	int    LatticeSize;
	double L0;
	bool   relax;
};
// *******************************************************************
// class lattice_1DTRI is comprised of 
// - state variable:
// - lattice variables;
// - Time variables
// - diffusion rate table
// *******************************************************************
class lattice_1DTRI
{
private:
	// Settings ----------------------------------------------------
	double	t_start;
	double	t_final;
	double	t_int;
	int		LatticeSize;
	int		M_ini;					// initial lattice layer number
	int		M_inc;					// increasement of lattice layer number
	double	temp_ini;
	double	dep_rate_ini;
	// INPUT -------------------------------------------------------
	double	Temperature;			// K
	double	Dep_rate;				// layer/s
	// OUTPUT ------------------------------------------------------
	struct lattice_state LState;
	// Internal variables ------------------------------------------
	int		Repeated;
	int		TopLyr;					// the top layer containing particles
	int		MLyr;					// the number of effective monolayers in the lattice
	double	dep_time;				// deposition time
									// when rec_idx equals record.record.rec_num, all records have been updated
									// rec_xxx[rec_idx-1] refers to the current state, which is consistent with the current lattice
	
	double	diff_ttl;				// total diffusion rate
	double	DIFF_RT[7];				// diffusion rate table: NN (number of nearest neighbors) to rate
	unsigned int	**pp_lattice;	// pointer of pointer to the lattice
	// Dimension = MLyr
	double	*p_diff;				// pointer to a diffusion rate table	
	// Dimension = LatticeSize
	int		*PREV_ST;
	int		*NEXT_ST;
	int		*hts;
	int		*hts_half;
	double  *db_hts;
	int		*occu;
	// Dimension = MODE*(2*LatticeSize)
	double	*phi_nx;
	double	*psi_nx;
	
	// WRITE TO FILES -----------------------------------------------
	// Average record
	struct Recorder{
		char*	FileName;
		int		rec_num;				// number of nodes in the state vectors
		int		rec_idx;				// index of the next record to be updated, ranges from 0 to record.record.rec_num-1
		struct	OneRecord *pRecords;	// Dimension = record.rec_num
	} record;

	// Detailed record
	char		FileName2[256];

	int			Detail;
	ofstream	fd_time;
	ofstream    fd_T;
	ofstream    fd_W;
	ofstream    fd_SOR;
	ofstream    fd_H;
	ofstream    fd_R2;
	ofstream    fd_M2;
	ofstream    fd_alpha;
	ofstream    fd_beta;

	// FUNCTIONS --------------------------------------------------
	void	init_lattice();
	void	init_diff_rate();
	void	init_p_diff();
	void	init_phi_psi();
	void	init_prev_next();
	void	initiate_stats();
	
	// Lattice operation
	void	fin_lattice();
	int		inc_lattice();
	void	save_lattice(int repeat, int Recorder_tIDX);
	int		test_lattice_cons();

	// Key KMC simulation functions
	void	get_relative_sites(int lyr, int st, int *nb_lyrs, int *nb_sts);
	void	get_relative_site(int lyr, int st, int nb, int *nb_lyr, int *nb_st);
	
	void	add_particle(int lyr, int st);
	void	remove_particle(int lyr, int st);
	int		add_relax(int lyr,int st);
	void	deposition();
	int		diffusion();
	double	OneMCTime(double r_p);
	
	// Output functions
	double	calc_SOR();
	void	calc_rough();
	void	update_records(double iCurTime);
	void	avg_stats();
	
public:
	// 2009.07.29 Gangshi: dep_time is protected and can be visited only by get_dep_time()
	// double	dep_time;

	// 2009.08.01 Gangshi: Add iDetail to regulate the output, 0: no detail files; 1: detail files
	lattice_1DTRI(int iLatticeSize, double iTimeEnd, double iTimeInt, double iTempIni = 850, double iDepRateIni = 0, int iM_ini=M_INI, int iM_inc = M_INC, char* iFileName = "lattice_1DTRI.dat", int iDetail = 0);
	~lattice_1DTRI();
	
	// Output routines
	void	get();
	double	get_dep_time();						// return current dep_time
	double  getTime();
	double  getH();
	double  getRho();
	double  getR2();
	double  getM2();
	void    getAlphas(double* oAlphas);
	void    getBetas(double*  oBetas);
	void    getSnapshot(const string filename);
	struct	lattice_state get_lattice_state();	// return lattice state at the latest recording point
												// if this state is use in the control design, it is a state feedback control
												// sampling delay may be added in the implementation by achieving previous data point
	void reset();			// must be called before running the first kMC simulation loop
	// 2009.07.30 Gangshi: introduce DepDuration as the duration of DepRate as a constant input to the kMC simulation
	void update(double iTemperature, double iDepRate, double iDepDuration);
	void update2(double time, double iW, double iT);

	void display_lattice();
	int  notfinish();
};

double	ranO(int *idum);
double	RandRand();
double	RandRandN1();
double	RandRandN0();

#endif
